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Abstract 

We carry out a systematic numerical stability analysis of ZND detonations of Ma- 
jda's model with Arrhenius-type ignition function, a simplified model for reacting flow, 
as heat release and activation energy are varied. Our purpose is, first, to answer a ques- 
tion of Majda whether oscillatory instabilities can occur for high activation energies as 
in the full reacting Euler equations, and, second, to test the efficiency of various ver- 
sions of a numerical eigenvalue-finding scheme suggested by Humphcrys and Zumbrun 
against the standard method of Lee and Stewart. Our results suggest that instabili- 
ties do not occur for Majda's model with Arrhenius-type ignition function, nor with 
a modified Arrhenius-type ignition function suggested by Lyng-Zumbrun, even in the 
high-activation energy limit. We find that the algorithm of Humpherys-Zumbrun is 
in the context of Majda's model 100-1,000 times faster than the one described in the 
classical work of Lee and Stewart and 1-10 times faster than an optimized version of 
the Lee-Stewart algorithm using an adaptive-mesh ODE solver. 



1 Introduction 

In this paper, we carry out a systematic numerical stability investigation for detonation 
solutions of Majda's model with Arrhenius-type ignition function in the high-activation 
energy limit, at the same time testing and comparing various different techniques for nu- 
merical stability analysis. Our results should have application also to the effective design 
of numerical methods for more complicated detonation models. 

Majda's model |Mlj . a simplified model for reacting compressible gas dynamics, is often 
used as a testing ground for theory and numerical methods designed for application to the 
more complicated reacting compressible Euler equations [BMR] IM2j . A question posed 
by Majda in [M2] is whether this simplified model is sufficient to capture the complicated 
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Hopf bifurcation/pulsating instability phenomena that occur for the full equations in the 
high-activation energy limitj^] To explore this theoretical question is one motivation for the 
present work. 

A second motivation comes from the numerics themselves. In their foundational paper 
|LSj on numerical stability analysis, in which they introduce the standard scheme now 
in use, Lee and Stewart describe the numerical determination of detonation stability as 
computationally intensive, and identify the development of more efficient schemes as an 
important problem in the physical detonation theory)^] 

We take advantage, therefore, of the simplified context of Majda's model, as a prov- 
ing ground for various different numerical schemes, in particular testing whether recently- 
developed techniques from the related problem of stability of viscous shock waves [HuZ2, 
IHLZ| |HLyZ[ IBHZl IZ3| can be imported in a useful way. 

Our main object is to test whether an alternative "Evans function- type" scheme pro- 
posed by Humpherys-Zumbrun [HuZlJ can outperform the standard scheme of Lee-Stewart, 
at the same time exploring optimal implementations for both. There is some reason to hope 
for improvement, since, as pointed out in [HuZlj. the Lee-Stewart shooting method can be 
reformulated as an adjoint Evans computation carried out in a backward direction, the di- 
rection from x = toward x = — oo in which eigenfunctions (normal modes) are required to 



decay. (See also Section 5.2 ) As pointed out in |Brt IBrZ| IHuZ21 IZ3j . there is a numerical 
advantage in shooting, rather, in the direction x = — oo toward x = that eigenfunctions 
are expected to grow, since error modes then decay exponentially relative to the mode being 
computed, with the advantage being proportional to the spectral gap between growing and 
decaying modes of the eigenvalue ODE. 

Other novelties of our investigation are the systematic use of adaptive-step mesh both 
in spatial and frequency variables, the development of higher-order high-frequency asymp- 
totics, the derivation of rigorous if sometimes conservative bounds on the maximum size of 
unstable eigenvalues, and the introduction of a hybrid and limiting schemes designed for 
the singular, square-wave limit. 



1.1 Results 

Our results are, first, that Majda's model does not appear to support instabilities of any 
kind for Arrhenius or modified Arrhenius ignition function, even in the high-activation 
energy limit. (For zero activation energy, see the analytical proof of |JYj .) 

Second, it does appear that the optimum version found for the Humpherys-Zumbrun 
scheme outperforms the optimum version found for the Lee-Stewart scheme, by a factor of 

1 "In particular, ... there is the possibility of Hopf bifurcation to pulsating reacting fronts with an 
associated exchange of stability..."- problem 1, p. 25, |M2) . 

2 "Finally, we point out that even though our scheme is direct and easy to implement, complete inves- 
tigation of the various regions of parameter space is computationally intensive. Any equivalent or more 
efficient numerical method for computation of detonation should be considered a valuable contribution and 
such approaches are needed to further explore the parameter regimes of instability."- closing paragraph, p. 

130, B3J. 
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1 — 10 depending on frequency, and on the average perhaps 2 — 4. However, this difference 
is dwarfed by the (scheme independent) one obtained by using an adaptive-step ODE solver 
in the x-evolution, which yields improvement over the fixed-mesh solver originally used in 
|LSj by a factor of 100 — 1000. Likewise, in our implementation, the difference between the 
implicit z-coordinatization of [Er2|, ILS| , in which the profile is explicitly solvable as a func- 
tion of z, and the x-coordinatization in which the problem presents itself, is improvement 
by a uniform factor of 6 (apparently due to the cost of evaluating the profile; see Remark 

Q| >. 

The message is that the choice of numerical scheme requires a bit of care. For, with 
the wrong choices, performance can degrade by a total factor of as much as 2400-24000! 
With all the right choices, on the other hand, the computation is for reasonable values of 
activation energy quite numerically well-conditioned, at least in the simplified context of 
Majda's model, with computation times on the order of that seen for a scalar or 2 x 2 
viscous conservation law, of a few seconds for an entire stability computation for a given 
set of model parameters. 

In the high-activation energy limit, as for any singular limit, the computational perfor- 
mance degrades. In such cases, we find it necessary to factor out as much growth/decay of 
the solution as we possibly can, with less effective schemes not even converging for reason- 
able precisions and computation times. 



1.2 Discussion and open problems 

As observed by Majda |M2j . the Majda model may be derived from the full reacting com- 
pressible Euler equations via weakly nonlinear geometric optics in certain limiting regimes. 
However, as noted in |Zlj . these can be seen to lie in the small heat-release/small-activation 
energy region for which detonations are known to be stable [Z4]. Thus, there is no phys- 
ical reason to expect that parallels should extend to the high-activation energy limit and 
associated pulsating instability phenomena. 

Nonetheless, the structural analogy between the equations persists, and so the observed 
results of universal stability are somewhat striking. It would be interesting to pursue further 
the difference between the two models, both at formal and rigorous levels. In particular, a 
very interesting open problem is the analytic verification of stability for general £, as done 
for £ = in [JY] , 

A novelty here as regards detonation literature is the derivation of rigorous bounds on the 
size of unstable eigenvalues. However, these in some cases degrade rapidly in the singular, 
high-activation energy limit, blowing up as as £ — > ooj^] Indeed, our convergence 

studies show these bounds to be much too conservative, suggesting a sharp bound rather 
of 0(£). The application of semiclassical limit/turning point techniques to obtain better 
bounds would be a very interesting technical question; see Remark |6.6| Another very 
interesting open problem would be to carry out a complete stability analysis in the limit 
as £ — y oo, a problem intermediate to the analysis of general £. This together with the 



3 In some interesting cases, including the square-wave limit, they do not; see Remark 6.7 
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numerical studies carried out here for bounded £ , would resolve the question of general £ 
by a combined numerical and analytical approach, similarly as was done for viscous shock 
waves in [HLZl |HLyZj iBHZllBLZ] . 



We have carried out here a systematic confirmation/examination in a simple context of 
a number of aspects of numerical detonation stability analysis, which we hope will serve 
as a useful reference for further developments. In particular, a very interesting direction 
for future investigation is to determine whether the gains in efficiency observed for Majda's 
model carry over to the full reacting compressible Euler equations. 

2 Preliminaries 

2.1 Equations and assumptions 

Consider the inviscid Majda model 

i 2 

(2.1) - ' V 2 



u t + [ — ) = qk(j)(u)z, 



z t = -k(f>(u)z, 

u, z, 4>,q, k £ M 1 , q > 0, k > 0, with Arrhenius type ignition function 



(2.2) <j>{u) 



' Ce -£/T(n) T>0j 
T < 0. 



Here T{u) is a relation approximating the temperature/velocity relation for a full ZND 
profile. We consider here the simple cases T(u) = u as proposed by Majda [BMRl IMT| IM2j 
and T(u) = 1 — (1 — 1.5) 2 , a downward quadratic relation qualitatively similar to that for 
the full ZND equations as proposed by Lyng-Zumbrun |LyZ2| . 



A strong detonation wave of (2.1) is a traveling- wave solution 



(2.3) (u, z)(x,t) = (u, z)(x — st), lim (u, z)(^) = (u±, z±) 

£— >±oo 

in the weak, or distributional, sense, smooth except at a single shock discontinuity, without 
loss of generality at x = 0, where u jumps from u* = u(Q~) to n(0 + ) as x crosses zero from 
left to right, with z_ = 0, z + = 1, > Uj > u+, and 

(2.4) u_ > s > u + , 

where Uj, defined as the minimum value of u for which T(ui) = 0, is the ignition temperature. 

That is, a strong detonation wave consists of a shock advancing to the right into a 
quiescent (i.e., nonreacting) constant state with reactant mass fraction z = 1, raising u 
above ignition level U{, followed by a smooth "reaction tail" in which combustion (reaction) 
occurs, in which z decays exponentially to value z_ = and u to 1 < u~ < u* as x — > — oo 
|EFlllMTip^FZTl|ESFZ2llJEWl lZIllZ2|. 
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2.2 Parametrization 



By the invariances of ( 2.1 ), we may take without loss of generality u* = 2, u+ = 0, s = k = 1, 



and C > arbitrary, leaving only the parameters q > and £ > 0. From (2.4) and the 
assumption that (u, z) converges as x — > dboo, we find |JY| IZlj that 

(2.5) (u, z){x) = (u+, z+) = (0, 1) for x > 



and also (see (2.7) below) u_ = 1 + y/1 — 2q, so that 

(2.6) 1 < u- < 2, < g < -. 

with activation energy £ varying in the infinite range < £ < +oo. 

2.3 Profiles 

The ZND profile equations may be written, adding q times the second equation to the first, 
as — s(u + qz)' + (%f) = an d — sz' + k(f)(u)z = 0. Integrating the first equation from -oo 
to x and solving the resulting quadratic, we obtain (see |JY} IZ1] for details) 

(2.7) u{x) = 1 + y/l-2q(l-z(x)). 
We then obtain z by solving the ODE 

(2.8) z' = k<j)(u(z))z, z(0) = 1. 

In the simplest case £ = 0, we have the explicit solution z = e kx . 

2.4 Square-wave structure and the high-activation energy limit 

In the high- activation energy limit £ — > 0, the profile ODE becomes singular, with vari- 
ation in z concentrated near the value u max for which the "temperature function" T(u) 
is maximized. For the Arrhenius case T(u) = u, this means concentration near the value 
Umax = u(0~) = 2 at the right endpoint x = of the profile, and results in a sharpened 
reaction spike near that point; see Figure [2] Here, following the standard normalization of 



|Er2, LS], we have chosen C in (2.2) so that z(—2) = 1/2, that is, the half-reaction occurs 
at a specified spatial point x = —2. For the modified Arrhenius case T(u) = 1 — (u — 1.5) 2 , 
rather, (j)(2) — > exponentially in £, and so the profile takes on a characteristic "square- 
wave" shape similarly as for the full reacting Euler equations, consisting of a long flat tail 
from x = — oo to an intermediate point xq for which u(xq) = u max = 1.5, a rapid change 
in the vicinity of xq, and a second long flat region from x = xq to x = 0; see Figure [TJ 
Here, we have have chosen a bit less carefully the normalization C = e 21 ^/ 40 in order to 
keep in frame the value of x at which z = 1/2. (Recall that change in C amounts to rescal- 
ing in x, so is not essential to our analysis.) For the full equations, square-wave structure 
and the high-activation energy limit are associated with Hopf bifurcation and transition to 
instability [EEO EH [LS] . 
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(a) 



(b) 



(c) 



Figure 1: Profile plots for fixed q value and £ = 1, 10, 20, 30, 40. We use a solid line for u and a dashed 
line for z plotted against x. Here 0(u) = Ce^ £/T(ll) , T(w) = 1 - (1.5 - u) 2 , and C = 10 21£/40 . In 
Figure (a), q = 0.15; (b), g = 0.3; (c), q = 0.45. 




Figure 2: Profile plots for fixed q value and £ = 1, 10, 20, 30. We use a solid line for u and a dashed line 
for z plotted against x. Here 0(u) = Ce~ £ ^ u , where C is chosen so that z(—2) = 1/2. In Figure 
(a), q = 0.15; (b), q = 0.3; (c), q = 0.45. 
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3 Linearized stability and the Evans— Lopatinski determinant 



We now briefly review the linearized stability theory of |ErH IJLW1 \Z1\ IHuZlj . Shifting 



to coordinates x = x — st moving with the background Neumann shock, write (2.1) as 
W t + F{W) X = R(W), where 

To investigate solutions in the vicinity of a discontinuous detonation profile, we postulate 
existence of a single shock discontinuity at location X(t), and reduce to a fixed-boundary 
problem by the change of variables x — > x — X(t). In these coordinates, the problem becomes 
W t + (F(W) - X'(t)W) x = R(W), x ^ 0, with jump condition X'(t)[W] - [F(W)} = 0, 
[h(x,t)] := h(0 + ,t) — h(0~,t) as usual denoting jump across the discontinuity at x = 0. 



3.1 Linearization/reduction to homogeneous form 

In moving coordinates, W is a standing detonation, hence (W°,X) = (W°,0) is a steady 
solution of the nonlinear equations. Linearizing about (W°,0), we obtain the linearized 
equations (W t -X'{t)(W°y(x)) + (AW) x = EW, with jump condition X'{t)[W°]-[AW] = 
at x = 0, where A := {d/dW)F, E := (d/dW)R. 

Reversing the original transformation to linear order, following [JLW] , by the change of 
variables W — > W — X(t)(W )' (x), and noting that x-differentiation of the steady profile 
equation F{W°) X = R(W°) gives (A{W )(W°Y(x)) x = E(W°)(W°)'(x), we obtain mod- 
ified, homogeneous interior equations Wt + (j4VF) x = EW together with a modified jump 
condition accounting for front dynamics of X'(t)[W°] - [A(W + X(t)(W )')} = 0. 

3.2 Evans— Lopatinski determinant 

Seeking normal mode solutions W(x,t) = e xt W(x), X(t) = e xt X, W bounded, of the lin- 
earized homogeneous equations, we are led to the generalized eigenvalue equations (AW)' = 
(-\I + E)W for x + 0, and ^(A^ ] - LA(T^ )']) - [AW] = 0, where "/" denotes al/dx, or, 
setting Z := AW, to 

(3.2) Z' = GZ, x^ 0, 



(3.3) 
with 

(3.4) 

and 

(3.5) 



X(\[W U ] - [A(W )'}) - [Z] = 



A 



u-1 



, E 



G = (E- XI)A~ 



qkd(p(u)z qk(p(u) 
—kdip(u)z —k(p(u) 



zMmi k(j) (u) + A 



u- 1 
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The Lopatinski determinant is then defined as 



(3.6) D ZND (A) :=det(Z- (A, 0), X[W] + R(0~))\ x=0 , 



where Z~(\,£) is a bounded exponentially decaying solution of (3.2), analytic in A and 
tangent as £ — > — oo to the subspace of exponentially decaying solutions of the limiting, 
constant-coefficient equations Z' = G-Z, and 

- ( !tS ) ■ ™ - ( "o 2 



More precisely, Z~ ~ e )J ~ x V-(X), where \i- is the positive real part eigenvalue of G_ and 
V- is an associated eigenvector chosen analytically in A [Br. IHuZ2l IZ3| . Evidently, there 
exists a normal mode solution with frequency A, 9iA > 0, if and only if Dznd(X) = 0. 



4 Numerical approximation 

To estimate Dznd numerically, we approximate Z~ at a large but finite negative value 
x = — M with the value Z~{—M) = e~^~ VL, and solve from — M to using a standard 
adaptive-step Runga-Kutta scheme. The vector V_(A) is evolved analytically by solution 
of Kato's ODE d\V = P'(X)V, where P (A) is the eigenprojection associated with see 
|BrZ| IZ21 IHuZ2l IZ3| . Numerical convergence and efficiency of general schemes of this type 
are discussed in |Br| IHuZ21 IZ3j . 

Remark 4.1. As discussed in \HuZl\ lZ^ , the single most important factor for efficiency 
of such computations is to use an adaptive-step rather than fixed-step ODE solver since 
traveling front and boundary-layer solutions inherently involve multiple scales. As seen 
here, failure to use adaptive steps can reduce efficiency by two or more orders of magnitude. 



Detection of roots. Zeros of the Evans-Lopatinski determinant can be found through 
individual A-evaluations by Newton's or other root-finding/following methods as in [LSJ. 
Alternatively, as in |Er2| \BiZ\ IHuZ2| IZ4j . they can be detected by a Nyquist diagram, or 
winding-number, computation, mapping a large semicircle S contained in the positive real 
part half-plane 3?A > via Dznd and taking the winding number to determine the number 
of zeros lying within S. We follow the latter approach here. We discretize S by adaptive A- 
steps, taking care that the relative change in DzND{Xj) is < 0.2 for each step, thus ensuring 
an accurate winding number count \Br\ IBrZj . 

z-coordinatization and renormalization. For purposes of numerical approximation, 
is advantageous to use the coordinatization of Erpenbeck, Lee-Stewart, and others, by z 
instead of x. This means replacing Z' = GZ by 

(4.1) Z = GZ, 
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where ' denotes d/dz and 



< 4 - 2 > : = ? = wr 

then integrating from z = to z = 1 instead of from x = — oo to x = 0. In practice, we 
integrate from z = e - fc< ?H- 00 ) A: f to z = 1, where x = — M is our usual starting point in 
a>coordinatized version, and initialize Z~ as usual as e~^ M V-. This avoids the need to 
solve the z equation z' = k<p(u)z numerically. 

A further (standard; see [Br|, IHuZ2] ) improvement is to renormalize Z = efo dividing 
out expected growth/decay, converting Y' = (G — fiI)Y to Y = HY, where 

(4.3) H{z) := ^ G — ill 



z' k<p{u(z)) 

Here, \i can be either /i_ or (better if computable analytically) the eigenvalue \i{x) of G(x). 

Remark 4.2. For an adaptive-step ODE solver, we found no mathematical difference be- 
tween x- and z- coordinates, as both required the same number of mesh points/functional 
evaluations for a given x (resp z) integration^ However, in practice, the change to z- 
coordinates gave an improvement of 6 times or more in speed due to the cost of the interpo- 
lation step used to evaluate the numerically pre- computed profile at the variable points needed 
for an adaptive step ODE. This could be improved somewhat by the use of a more efficient in- 
terpolation scheme; however, z-coordinates are always preferable for an adaptive-step ODE 
solver. As discussed in JZ31 \HuZlf . renormalization typically improves speed for a single 
X-evaluation by a factor of 2 or more. This is magnified in the high- activation energy limit, 
for which growth/decay rates become extreme, and the unrenormalized computation often 
cannot even be carried out, leaving machine scale and returning NaN errors. For winding 
number computations, it is still more important to divide out expected growth/decay, which 
otherwise introduces additional winding, decreasing the A stepsize and greatly increasing 
computational cost. 



Reduction by A. By translation-invariance, Dznd has always a root at A = 0. Like- 
wise, there is an extra factor A as A — > oo induced by the form of factor (A[l^°] + R) in the 
defining determinant beyond what is induced by growth/decay of the ODE solution Z~; see 
also the more detailed discussion of Section [6j For both these reasons, it is advantageous in 
performing winding-number computations to work with a reduced Evans-Lopatinski deter- 
minant Dznd(^)/ 'A, effectively removing a single zero, hence one circuit about the origin, 
and thereby greatly reducing the number of A-points required for the computation (typically 
by factor 4 or so); see Figure [3j All winding number computations are done, therefore, with 
respect to the reduced determinant, throughout the paper. 

4 Additional experiments, not recorded here. 
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Figure 3: For £ = 10, q = 0.3, G = io B / 2 + B / 40 , and R = 4, which we get from the high-frequency 
convergence study, it takes 2.7524 seconds and 127 mesh points to obtain 0.2 tolerance in the A 
contour for the unreduced Evans-Lopatinski determinant Dznd- If we divide by A, it takes 0.6858 
seconds and 28 mesh points. Open circles correspond to the reduced Evans function for which we 
divide by A. 

5 Alternative formulations 

5.1 The adjoint Evans— Lopatinski determinant 

We may define an alternative stability function (sec [HuZ2, CJLW, Z4J) as 

(5.1) DzndW = (Z, (X[W°] + R))\ x=0 , 

where Z denotes the unique (up to scalar multiplier) decaying solution of 

(— qk d(f>(u)z+\* kd<f)(u)z \ 

Dznd and Dznd differ by a nonvanishing analytic factor j3, by duality relations Z = -yr 
and det(u, w) = ■ w, where W = ce-^o ( G ) is an appropriately normalized Wronskian 



of (3.2) and D Z nd = (3D ZND with P = W. 



Remark 5.1. The method of reduction to homogeneous form JJLW\/ is equivalent to solving 
the inhomogeneous equation 

(5.3) Z' = GZ + XXW'(x) 

by linear superposition, using the fact that a particular solution is given by W p = XW'(x), 
hence W = W~ + XW'(x), or Z = Z~ + XAW'(x) = Z~ - XR{x), and substituting into 
the jump relation \X\W®\ + [Z] = 0. 



As discussed further in jHuZlj . (5.2) may be viewed as a streamlined version of the 
original method of Erpenbeck [Er2j; we shall refer to this scheme as the homogeneous 
Erpenbeck method. When renormalized by /2_ (resp. fi(x)), we will refer to it as the adjoint 
H (resp. adjoint /J,(x)) method. 
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5.2 The method of Lee and Stewart 



The method of Lee and Stewart consists, rather, of solving the inhomogeneous equation 



(5.3) with X = 1 from x = initialized with Z(0) := A[W]. They then take the inner 
product of V~ against Z(—M), with M chosen sufficiently largej^] By duality, this can be 
seen to be exactly 

(5-4) D LsW = Jf^, 

where /U2(— oo) is the negative eigenvalue of G-. 



Expanding W' = (dW /dz)z' = (dW/dz)k<pz, we may rewrite (5.3) in z-coordinates as 

/ du/dz\ 



(5.5) Z=^-Z + \( d ^ d2 



where du/dz is determined from the profile solution (2.7). Equivalently 
(5.6) Dznd(X) = e-^^ M D LS (X) = (V~, Z c (—oo)), 

where 

(5.7) 

Homogeneous version. The Lee-Stewart method may equally well be implemented 
via the homogeneous equations Z' = GZ, Z(0) = X[W] + R. (Here, we are using Z — Z — > 
as x — > — oo.) This appears to be slightly faster but essentially equivalent in practice. 

5.3 The Polar/Drury method 

An alternative to direct renormalization is the polar method (or continuous orthogonaliza- 
tion method of Drury) [HuZ2], which, in z coordinates, appears as 

(G - YY*G)Y (G - YY*G)Y 



(5.8) Y 



k(j){u{z))z 



where D po i ar (\) := det(Y, A[W" ] + R)\ z =i- This differs from Dznd by a nonvanishing 
continuous factor r(A), so has the same winding number and roots. It has the advantage 
that \ Y\ = 1, so growth/decay is completely factored out, while good numerical conditioning 
is maintained so long as there remains at least a neutral spectral gap between /x_ and the 
other eigenvalue of G- |Z3j , as holds in this case on a neighborhood of 5RA > (see Section 
rcjj). An adjoint version can be computed similarly; we refer to this scheme as the polar 
adjoint method. The factor r(A) may be computed by integrating logr = fc( ^ M ^-^T from 
z = to z = 1, with r(0) = 1, and used to recover Dznd |HuZ2]. We refer to the latter 
scheme as the polar radial, and the analogous dual as the polar adjoint radial method. 



5 More precisely, since they work like Erpenbeck in z and not x coordinates, against Z\ z=e , e > 
sufficiently small. This is phrased in terms of a progress variable A := 1 — z; what they call a is our A. 
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5.4 Hybrid and limiting methods 



In the high-activation energy limit £ — > with modified Arrhenius-type ignition function, 
u approaches the well-known "square-wave" profile familiar from the full ZND equations 
\FD\ lErlj . which consists approximately of three constant zones separated by a reaction 
jump following at some distance the Neumann shock. This singular configuration is difficult 
for the standard numerical Evans-Lopatinski function to resolve, since the integration passes 
through an unexpected second, different, constant region for which it is not tuned. Indeed, 
following the numerical prescription of [Z3J, we should for best numerical conditioning 
rather perform an Evans function type computation integrating from — oo forward toward 
the reaction jump and from backwards toward the reaction jump, taking a determinant 
in between. 

This can conveniently be accomplished by a hybrid method intermediate to the adjoint 
and Lee-Stewart methods, in which we compute from z = to z = 1/2 by the adjoint 
method, in some optimized version, and from z = 1 to z = 1/2 via the Lee-Stewart method, 
again in optimized version, taking the inner product at z = 1/2. By duality, this gives the 
same result as the adjoint method with same optimized formj^] This should be done only 
when (ft (2) « ma,x(ft, so that such singular square-wave structure occurs. When this is not 
the case, the adjoint method is expected to be preferred. However, the hybrid method is 
safe, in that its numerical conditioning should always lie somewhere between that of the 
adjoint and Lee-Stewart methods, and captures some of the gain of the adjoint method 
even when square-wave structure does not occur. 

Limiting method. When (ft (2) « max (ft, and square-wave structure is truly in effect, 
we may take a limit before computing the hybrid determinant and simply project out decay- 
ing modes. That is, we may define a projective, limiting, method, in which we substitute 
for initial data X[W] + R the data P(A)(A[W] + R), where -P(A) is the eigenprojection of 
G\ z= i associated with the negative real part eigenvalue n%. For related analytical results in 
a similar situation, see |Z5j . We will refer to this version as the Evans function method. 



6 High frequency asymptotics 



Making the change of coordinates x — > x := |A|x, A — > A := A/|A|, we convert (5.2) to the 
approximately diagonal system 

(6.1) Z = A(ex)Z + eB(ex)Z, 

where |A| = 1, 5ftA > 0, e := |A|~ , ' denotes d/dx, and 

/ A* —egk d<j>{u)z n \ / fl kd<p{u)z \ 

(6-2) A =\ V , „_J. B 



u-1 



-A* - ek(ft(u) I ' \qk(p{u) 



6 In practice, we solve Z' = (G — fj,2{x))Z, where ^ is the negative eigenvallue of G, and the adjoint 
version Z' = -{G - fi 2 (x))*Z. 
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Making the further change of coordinates Z = TY, 



kd(j){u)z qk(f)(u)(u - 1) 

(6.3) a = — — = o = 



— itA* — e(k(f>(u) + qk d<p(u)z) uX* + e(k(f>(u) + qk d(j)(u)z) 

we obtain 

(6.4) Y = T~ X {A + eB)TY - T~ lr TY = A X Y + e 2 B 1 Y + e 3 S 2 y, 



where 

(6.5) A 1 



( Oil 


) 


■{ 




OL2) 





A* —eqk dif>(u)z ^ 
u-1 U 



-A* - ek(p(u) 



(6.6) 



-aqk^u) + - a&(A* + k<t>{u)) -a'qk^u) + ^ 

bx i_ b 2 kd<j>(u)z n r,hrh(T;\ _ bk d 4>(u)z _ „ u \* -qk d<t>(u)z 



Bo 



ab x q 



1— e 2 a6 

1— e^ao , 



6.1 Estimation of Z 

Change now to y-coordinates, where dx/dy = 4>(u(x), i.e., dz/dy = ekz, or z(y) = e ke ^. 
Defining Y(y) = e^o a ^i z ) dz V(y), and simplifying, we obtain as usual the Duhamel repre- 
sentation: 

(6.7) V{y) = TV{y) := VI + f e ^ ^- &lI ^ d ' s {e 2 B x + e 3 B 2 )V(z)dz, 

J —CO 

where A = A/<f>, aj = otj/cj), B 3 = Bj/<j>, and \B x {y)\ < Cie~ ke \y\, \B 2 (y)\ < C 2 e~ k£ M for 
some Cj > 0, and V- = V^_(A) = c(A)(l,0) T is the asymptotic limit at x — > -co, c / 0. 
Here, by the Mean Value Theorem, Cj < sup \dBj/dz\. From a 2 — a\ = ~" A ~ gfc ^(yH £ g fcrf( ft(") z ) _ 
we have by positivity of ek(f>(u-), 3?Au, and u — 1, the bound 

(6.8) ^//(Ai-diX*)^! < e ft* t <> &2 - &1 ')W dS < e 7 ^ S d2 < e rz{ y ] < e 7 , 

(6.9) ,:=,»p(Pf), 

x \ U — 1 / 

where denotes the positive part of /. For 3?A > max ^fcM ) we have simply 

(6.10) i e ;/(^i-ai)(*)«B| < e ;/»(fi2-ai)(i)ds < L 
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Remark 6.1. For the Arrhenius-type ignition function <fi(u) = Ce^ , we have 

(6.11) dcf)/(t) = £/u 2 > 0, hence 7 = q£/u1 < q£. 

For the modified Arrhenius ignition function (p(u) = Ce T (") ; where T{u) = 1 — {u — 1.5) 2 , 
we have 

(6.12) d(j)/(j) = -2{£/T 2 ){u - 1.5), hence 7 < (4/3)g£ for q > .375 
and 

(6.13) 7 = for q < .375, in which case u > u_ > 1.5 and d<f>/<f> < 0. 
A sufficient condition for KA > max iMi^zhi.^ i s 

(6.14) 9iA > /c7max(/» ~ kSmaxcf). 

Lemma 6.2. For e < 2e - l{c k +eC2) or e < 2(Cl * eC2) and KA > max qkd4, ^~ k<t> , T is a 
contraction on L°°( — 00, 0], and Z~\% = q is proportional to (l,^) 7 where lo<i := 
and \u\ < £ -f(C l + sC 2 ). 



Proof. From (6.7) and the stated bounds, the Lipshitz norm of T is bounded by 

-2/ 



^ {C \t £C2) e~ k ^dz < ee*(C x + eC 2 ) f 

J-00 n z J-c 



( 6 - 15 ) <ee\C 1 + eC 2 ) / e~ ke ^dz 



e 

00 
y 



-ke\z\ d ~ 



OC 



< ^(Ci + eC 2 )e-*l*l, 

hence it is contractive for ^-{C\ + eC2) < 1/2, or e < 2e^(c , H-eC7 2 ) ' w ^ n relative error at 
rc = from the diagonal solution (i.e., with e set to zero) given by ^-{C\ + eC 2 )- 

That is, Y(0) is proportional to (1,uj) t , with < tt(Ci + eC-i). Transforming back 

by Z~ = lY, with T = \ , we thus find that 

-1— (A t) C) - 

is proportional to (l,w 2 ) T , where w 2 := ^^/q v, ■ □ 

Remark 6.3. Lemma \ 6.S\ may be recognized as a variable- coefficient analog of the gap 
lemma of JGZf . and a quantitative version of the abstract Lemma A.l established for the 
full ZND system in [ZJ^]. Though we did not state it, the argument implies also that Z~ 
converges exponentially in relative error as x — > —00 to the solution of the diagonal system 
Z' = A\Z , A\ as in (|6.5[). A similar argument applies for the related polar method t r HuZ2j. 
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Corollary 6.4. D ZND {\) / for < KA and 

|A| > max{— (C x + eC 2 ), -^(q + |w 2 |)}, w 2 



2 l/J ' ' 1 + ea(0)w ' 

or 

qkd(pz - k(j) k(p(2) eb(0) + u 

MX > max = ana |A| > (g + |w 2 |), ^2 ~ 



u 1 1 ~ 2 ^ 1 1 + ea(0)w 

Proo/. Computing A[W] + R(W(Q-)) = (-2A + qk<p{2), -k<p{2)) T , we have that D ZND {\) 
is proportional by a nonvanishing analytic factor to 

(Z-(A,0), (-2A + ?fc0(2), -^(2)) T )U. =0 = (1, w 2 ) T • (-2A + g fc0(2), -A:0(2)) T , 

and thus to (l,w 2 ) T • (— 2A + eqk</)(2), — ek(p(2)) T = — 2A + £&0(2)((j[ — w 2 ), which is nonva- 
nishing provided ek(f>(2)(q + |w 2 |) < 2, or e < fc0(2)( g +M) - D 

Remark 6.5. When \X\ > ^(C a + eC 2 ), |w| < 1, so that \ u 2 \ is roughly 1 as well, and 
k <t>( 2 ) ( i I i\ <r 3fe0(2) 

Remark 6.6. Corollary \6.J\ is efficient for X near the real axis. However, for the Arrhenius- 
type ignition function T(u) = u it can be quite inefficient away from the real axis, where it 
does not make sufficient use of spectral separation of modes as opposed to spectral gap, hence 
requires the unusable bound \X\ > e 7 ~ e £ as £ — > oo. To obtain practical (guaranteed) 
bounds by use of more modern turning-point theory is an important direction for future 
investigation^ 

Remark 6.7. For the modified Arrhenius ignition function (quadratic version), our esti- 
mates are efficient for q < .375, where 7 = 0: in particular for the "difficult" case q = 0.3 



we have much investigated; see Remark 6.1 



Remark 6.8. For q bounded from 1/2, k = 1 and £ = 0, C = 1, we have 7 = 0, a = 0, and 
b, b x roughly 1, giving C 2 ~ C\ ~ 1. Thus, e < 1 is sufficient in this case for contraction, 
or \M ^ 1- Our bounds blow up in the CJ limit q — > 1/2, for which (u — 1) — >• as x — > —00. 
(However, note that 6—^0 as q —> q c j, partly compensating for badness of this limit.) This 
case would be interesting for further investigation. 

6.2 Limiting behavior 



By arguments like those of Section 6.1 and especially Remark 6.3 we find that 
Z"(0) 



' J-^{fc{v)-p2(~°°)dv\ ( e kf° 00 (<f>(u(y))-<t,(u{-oo)))dy S 



/ V 



as e — > 0, where /3 2 := — a\ = A + ek(f>(u), a 2 as in (6.5), yielding the following asymptotic 
behavior. 



7 For the Arrhenius-type ignition function, actual instabilities for £ — ¥ 00 are expected to appear (if 
indeed they do) only in the region |A| < max qkd(f>z ~ _Emax</> where spectral separation |qx — » 1 
fails. 
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Proposition 6.9. For 5ftA > 0, 

(6.16) D ZND (X) = CX(l + 0(\- 1 )), D ZND {\) = C l \e c * x {l + 0{\~ 1 )). 

as |A| -> oo , where C := 2e fc ^-- Ci, and C 2 are real constants. 

Remark 6.10. Higher order approximants D Z ndW = e ClX+Co+c - lX ~ 1 ' ) A(l + 0(A~ 2 )), 
etc., may be obtained by further diagonalizations as discussed in \Z4\ IM"aZ3j . 

In particular, D ~ 2A for £ = 0. These relations can be used to determine a maximum 
radius through a convergence study, which, as pointed out in [HLZ] HLyZ| , in practice 



typically gives better bounds than those obtained by rigorous tracking/conjugation bounds. 

7 Numerical experiments 
7.1 Initialization at — oo 



In our numerical experiments we compute the Lopatinski determinant given in (3.6). Fol- 
lowing \L2\ IZ3| , we use the ODE 

(7.1) S' = P'S 

of Kato to explicitly compute an analytically varying initializing eigenvector V~ for the 
Evans function, where ' denotes d/dX, 5(A) is the desired eigenvector of G_oo, and P(X) is 
the associated eigenprojection. Preserving analyticity in this manner, we are able to employ 
the argument principle thus determining the number of zeros of the Evans function inside 
a contour by computing the winding number. At x = — oo, 

where a = = ^/\-2q ' ^ = ~lk4>- } and c = k(j)—. Then right and left eigenvectors of 
G-oo corresponding to the eigenvalue c + A, which satisfies 5R(c + A) > for 5RA > 0, are 

/ -b \ r 
respectively r+ = i)*—c I / + = (o l). We then have the projection P = 



1 



-b 



it's derivative (with respect to A) P', P = ^ (a-i)\-c I ; p/ 
the ODE 5' = P'S, we find that 

/ -6 \ 




(7.3) S(\) = a 2 \^ a - 1 l x - c J+a l ^ ) 

where we take a\ = and Q2 = 1. The corresponding initializing vector for the adjoint 

- / 1 \ 
method is V- = & 

\ ( a -l)A*-c / 
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7.2 Computation of Z /evaluation of D 

Following the general approach of [HuZ2| iHLZf IBHZ] . the ODE calculat ions for individual 
A are carried out using MATLAB's ode45 routine, which is the adaptive 4th-order Runge- 
Kutta-Fehlberg method (RKF45). This method is known to have excellent accuracy with 
automatic error control. Our standard error tolerance setting is AbsTol = le-6 and RelTol 
= le-8 unless otherwise mentioned, and the value of approximate minus spatial infinity — M 
is determined experimentally by the requirement that the absolute error 



\(u,z)(-M)-(u,z) 



< TOL 



be within a prescribed tolerance, say TOL = 10 3 , where (u, z) and (u, .z)-oo are the pro- 
file and limiting endstates of Section 2.3 In the rescaled z coordinates where we do not 



even need the profile to compute, we perform a convergence study of the Evans-Lopatinski 
determinant for A = 2i to determine numerical infinity, requiring relative error between 
output be less than 10~ 3 between successive computations for numerical infinity decreasing 
as negative powers of 2. For a theoretical convergence analysis in the simple case consid- 
ered here, see |Brj : for a more general treatment, see |Z3| . For the results of a numerical 
convergence study, see Fig. [4] below. 

Fixed-mesh Lee-Stewart implementation. For comparison purpose, we carry out 
also experiments using a fixed- mesh grid as prescribed in |LSj . taking e = 1/N, and using a 
uniform mesh for z 6 [0, 1] of mesh points Zj = j/N, then choosing N large enough to get 
a prescribed level of convergence, as determined by numerical convergence study. 




(b) 




Figure 4: a) Convergence of the profile, b) Convergence of Evans approximations (regular method with fi(x) 
scaled out) on a contour with R = 1, for <f> = f , q = 0.3, k = 1 and approximate spatial infinity 
M = 0.9, 0.5, 10"\ 10" 2 , 10" 3 . Run times were respectively 1.08, 0.40, 0.76, 1.06, 1.33 seconds. 



7.3 Verification: comparison with exact solution 

For comparison, we test our code against an exact solution found in [JY] for 



1 of 



(7.4) 



Dznd(X) = (2A + (2 - q- qM)) (j^) 
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where 



(7.5) * 



_ r° 

e ^ 



-2,(1- 



a/1 - 2g(l - e k v) 



dy, 



Pit) 



a/1 - 2q(l - e fc ?) 



This formula results from a choice of Z~ asymptotic to e^ A+fe ^(*, 1) T as cc — )■= oo in 
agreement with (7.3), hence should agree with the results of our code. As seen in Figure [5j 
the agreement of code with exact formulat is excellent. 




Figure 5: Comparison in the complex plane to the exact solution of JY evaluated on a semicircle of radius 
R = 10. a) full size, b) zoom top right. Open circles correspond to the exact solution and the line 
with closed points to the fi(x) method. Here £ — so that tj> = 1, k = 1, and q — 0.3. 



7.4 Winding number computations 

To check stability, finally, we determine the number of roots within a semicircle S = 
dB(0,R) n {!RA > 0} of large radius R by computing the winding number of the image 
curve D(S) as 5 is traversed counterclockwise, which, by the Principle of the Argument, is 
equal to the number of roots within. We compute this winding number by varying values 
of A around S along 20 points of the contour, with mesh size taken quadratic in modulus 
to concentrate sample points near the origin where angles change more quickly, and sum- 
ming the resulting changes in arg(D(A)), using ^slogD(X) = argD(A)(mod27r), available in 
MATLAB by direct function calls. As a check on winding number accuracy, we test a pos- 
teriori that the change in argument of D for each step is less than 0.2, and add mesh points, 
as necessary to achieve this. Recall, by Rouche's Theorem, that accuracy is preserved so 
long as the argument varies by less than ir along each mesh interval. 

Computations were carried out within the MATLAB-based STABLAB code developed 
by J. Humpherys with help of the authors. Using MATLAB's parallel computing toolbox 
on an 8-core Power Macintosh workstation, we were able to achieve a speedup of over 600%, 
similarly as in the previous numerical studies |BHZl IBLZ| IBLeZ] . 

Determination of R. The radius R is chosen so large that there exist no unstable roots 
outside B(0,R), ensuring that the roots found by our winding number computation are the 
only possible ones on the nonstable half-plane KA > 0. This could be done analytically as 
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described in Section [6] (and carried out, for example, in IH.vX|). Here, following instead 

to match 



6.9 Remark 6.10 



|BLeZl IBHZj , we use the asymptotics described in Proposition 
the reduced Lopatinski determinant D(X)/X to a first-order or higher-order approximant 
e CA Qr e CiA+Co+C-iA carr yi n g out a convergence study to determine when D(X)/X has 
sufficiently converged. 

This was accomplished by, first, requiring that the relative error between D(X)/X and a 
best-fit value of Ce Cl ^ x+C2 ^ x be less than or equal to 0.2, and, second, that doubling the 
radius results in reduction of the error by a factor of approximately two, in accordance with 
the linear error dependence predicted by analytical theory. 

Typical convergence studies are illustrated in Figure [6] and Tables [T] and [2} A comparison 
of the Evans function vs. reduced Evans function is given in Figure [3j 



(a) 




(d) 




(<0 




(c) 




(/) 




Figure 6: Here we use solid dots to plot the Lopatinski determinant evaluated on a quarter arc, and we use 
open circles to graph the approximating function, C\ exp(C2/A). For the fi{x) method, we have 
Figures a-c corresponding respectively to R — 2, R — 4, and R — 2048. For the polar method, we 
have Figures d-f corresponding respectively to R — 4, R = 16, R = 64. Parameters are £ = 10, 
q = 0.3, <t>{u) = Ce-*™, C = 10 21£ / 40 . 



Results. We computed, using the Evans-Lopatinski determinant (our basic algorithm) 
rescaled by fJ>(x), a batch job for the parameter values, {£,q} = {0.01 : 0.01 : 0.49} x{0 : 0.1 : 
5, 5.2 : 0.2 : 10, 12, 15, 20, 30, 40}, with the exception of a few numerically challenging 
parameters {E, q} = {20} x {0.49}, {30} x {0.48, 0.49}, {40} x {0.47, 0.48, 0.49}, with C = 
lO 21 ^/ 40 , for ignition function jb(u) = Ce~ £ ^ T ^ where T(u) = 1 — (u— 1.5) 2 . Computational 
statistics are given in Table For ignition function <f>(u) = e s / 2 e — £ / u we computed the 
Evans-Lopatinski determinant for {£, q} = {0 : 0.1 : 5, 5.2 : 0.2 : 10, 12, 15} x {0.01 : 0.01 : 
0.37, 0.375, 0.38 : 0.01 : 0.49} U {20} x {0.01 : 0.01 : 0.37, 0.375, 0.38 : 0.01 : 0.47} U {25} x 



8 For winding number computations, we find that it is important to divide out as much behavior as we 
can, to avoid excess winding. Dividing by A is crucial, CieC2A still better. Eventually, there is a break-even 
point in complexity (and coefficient size of remainder) vs. power of 1/|A|, here occurring at first order. 
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Radius 


Relative Error 


KO 


Kl 


K2 


2 


0.00387783 


0.693873 


-0.157136 


0.0202546 


4 


0.000656158 


0.693278 


-0.152452 


0.0110315 


8 


8.88927c-05 


0.693167 


-0.150713 


0.00427183 


16 


9.83384c-06 


0.69315 


-0.150192 


0.000241910 


32 


6.99413c-07 


0.693148 


-0.15005 


-0.00195688 


64 


1.26154c-07 


0.693147 


-0.150013 


-0.003106 


128 


5.78215c-08 


0.693147 


-0.150003 


-0.00369357 



Table 1: For £ = 10, q = 0.3, C = e £ / 2 , and 4>{u) = e ~ £ / u we examine the high frequency convergence of the adjoint 
fi(x) method with it's best fit with e K0+Kl/X+K 2/A 2 ^ 



Radius 


Relative Error 


K0 


Kl 


K2 


2 


0.0436173 


0.72503 


-0.501318 


0.0641958 


4 


0.0294181 


0.706568 


-0.347815 


-0.254427 


8 


0.0100462 


0.696782 


-0.190304 


-0.88822 


16 


0.00172202 


0.693867 


-0.0981737 


-1.61602 


32 


0.000146306 


0.693263 


-0.0604492 


-2.20532 


64 


6.21728e-05 


0.693164 


-0.0480715 


-2.58952 


128 


2.54483o-05 


0.693149 


-0.0444984 


-2.81055 


256 


7.81771e-06 


0.693147 


-0.0435363 


-2.92936 


512 


2.14819e-06 


0.693147 


-0.0432867 


-2.99093 


1024 


5.61212c-07 


0.693147 


-0.0432227 


-3.02249 



Table 2: For E = 10, q = 0.3, <j>(u) = Ce- £ ' T{ - U \ T(u) = 1 - (1.5 -u) 2 , C = 10 21£ / 40 we examine the high frequency 
convergence of the adjoint method scaled by fJ,(x) to it's best fit with e K 0+Kl/ X+K2/ X 



{0.01 : 0.01 : 0.37, 0.375, 0.38 : 0.01 : 0.45} U {30} x {0.01 : 0.01 : 0.37, 0.375, 0.38 : 0.01 : 
0.40}. Computational statistics are given in Table |4j All computations yielded winding 
number zero consistent with stability. 



q 


E=0 


E=20 


E=30 


E=40 


0.01 


(le-4,4,10,0.0086,1) 


(le-4,40,10,0.018,lel) 


(le-4,60,10,0.023,lel) 


(le-4,80,10,0.026,2el) 


0.1 


(le-4,4,10,0.1,1) 


(le-4,40,21,0.2,3) 


(le-4,60,33,0.19,3) 


(le-4,80,40,0.2,4) 


0.2 


(le-4,4,11,0.15,1) 


(le-4,40,39,0.19,3) 


(le-4,60,54,0.2,4) 


(le-4,80,74,0.2,7) 


0.3 


(le-4,4,12,0.17,2) 


(le-4,40,49,0.19,3) 


(le-4,60,75,0.19,6) 


(le-4,256,89,0.19,2el) 


0.4 


(le-4,4,17,0.16,2) 


(le-4,40,62,0.2,4) 


(le-4,60,84,0.2,8) 


(le-4,512,107,0.2,4el) 



Table 3: For cj> = Ce~ £/T{u) where T(u) = 1 - (u - 1.5) 2 and C = e 21£40 we record statistics for our Evans- 
Lopatinski determinant computations. The data represented are (L,R, error, time) where L is the 
numerical value of infinity in the spatial z domain, R is the radius of the domain contour, error is 
the maximum relative error between A contour output, and time is the time it took to compute. 



8 Performance comparisons 

In the remainder of the paper, we collect a number of data comparing performance of the 
various methods. 

In Figure [7] we demonstrate agreement between methods that should yield the same 
output providing a verification of the correctness of our code. This also provides a visual 
example of which methods require a tighter A mesh in order to obtain relative error within 
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E=0 


E=10 


E=20 


E=30 


q=0.01 


(0.1,4,0-004,0.4) 


(0.1,20,0.0063,0.8) 


(0.1,40,0.0075,1) 


(0.1,60,0.0081,2) 


q=0.1 


(0.1,4,0.045,0.7) 


(0.1,20,0.073,2) 


(0.1,40,0.088,4) 


(0.1,60,0.096,6) 


q=0.2 


(0.1,4,0.1,1) 


(0.1,20,0.18,3) 


(0.06,40,0.11, le+01) 


(0.06,60,0.13,2e+01) 


q=0.3 


(0.1,4,0.18,1) 


(0.06,20,0.17,7) 


(0.06,40,0.18,2e+01) 


(0.03,60,0.18,7e+01) 


q=0.4 


(0.06,4,0.16,2) 


(0.06,20,0.15,le+01) 


(0.03,40,0.17,8e+01) 


(0.03,60,0.19,4e+02) 



Table 4: For cj> = e £ / 2 e~ £ / u we record statistics for our Evans-Lopatinski determinant computations. Here 
we use the polar adjoint method with ode solver tolerance set at 10 -12 . We set the radius of the 
domain contour to be the max of 2£ and the radius obtained by the best curve fit of the Evans- 
Lopatinski determinant with e E o+ E i/ x + E 2/^ The data represented are (L,R,error,time) where 
L is the numerical value of infinity in the spatial z domain, R is the radius of the domain contour, 
error is the maximum relative error between A contour output, and time is the time it took to 
compute. 



tolerance. 

In Table [6] we demonstrate computational time for the various methods for a fixed 
A mesh. In this experiment, we simply set tolerance at TOL = 10~ 12 for relative and 
absolute truncation error in the adaptive RK45 ODE solver and integrated around a fixed 
radius 10 semicircular contour in the complex A-plane of the type used in our winding 
number computations for stability, with a small semicircle of radius 10 removed around 
the origin. The time taken to complete this computation gives a rough "average" measure 
of performance over different A regimes. However, it is a bit conservative, as it measures 
truncation and not convergence error, so does not reflect the expected better numerical 
conditioning of the Humpherys-Zumbrun vs. other schemes. For this reason, we have 
performed a series of more careful experiments at individual A- values, specifying convergence 
rather than truncation error, described below. Note finally that the standard fixed-mesh 
Lee-Stewart method described in |LSj is not included in the comparisons of Table |6j since, 
not being adaptive, it does not allow specification of truncation error in this simple way. 
The important comparison to this scheme is carried out in the convergence error-based 
study below. 

For the convergence error-based study, we evaluate the various methods at the represen- 
tative values A = 1, i, 10, lOi and determine the relative tolerance with which the methods 
must be solved in order for the output to converge to a specified tolerance. For the meth- 
ods using an adaptive mesh in z, we specify absolute tolerance to be 10~ 15 and increase 
relative tolerance by powers of 10 starting at 10 _1 until consecutive output is less than 
10~ 6 in relative error. We use MATLAB's ode45 function which employs the fourth order 
Runge-Kutta-Fehlberg method. By specifying absolute tolerance to be 10~ 15 , we ensure 
that relative rather than absolute tolerance determines the mesh setting^] 

For the fixed z mesh method of Lee and Stewart, we determine a number of mesh points 
N such that relative error between the fixed mesh inhomogeneous method and the adaptive 
inhomogeneous method, solved with error tolerance requirements of 10 -10 , is less than 10~ 6 

9 The ode45 algorithm requires that relative tolerance or absolute tolerance be met. 
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Figure 7: We compare different methods of integration for £ = 10, q = 0.3, 4>{u) = Ce~ £ / Tl ~ u \ T(u) = 
1 — (1.5 — u) 2 , and C = 10 21£ ^ 40 , R = 10, with 65 fixed domain mesh points. These comparisons 
serve as code verification, showing agreement between methods that should give the same output, 
and contrast the number of points needed by the various methods to obtain output with relative 
error within tolerance, (a) polar adjoint radial method, (b) adjoint fi method, (c) adjoint n(x) 
method, (d) polar adjoint method, (e) polar method, (f) (i(x) method, (g) polar radial method, 
(h) fx method, (i) no fi method, (j) homogeneous Erpenbeck method (k) the homogeneous method 
of Lee and Stewart 
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but greater than 9 * 10 7 . For the fixed z mesh method of Lee and Stewart, the reported 
tolerance is the relative error between the ODE output of these two methods. 

In Tables [7}j24 we report the computational statistics for the n(x) method, the adjoint 
n(x) method, the hybrid method, the method of Lee and Stewart (adaptive mesh in z), 
the method of Lee and Stewart (fixed mesh in z), the polar adjoint method, and the polar 
adjoint radial method, for both ignition functions (f>(u) = eT^l u and (f>(u) = e~^l T{ ^ where 
T(m) = 1 — (u — 1.5) 2 . For the first ignition function we record data for S = 1, 10, 20, and 
q = 0.1, 0.3, 0.47, and for the second ignition function we give statistics for £ = 0, 10, 40 
and q = 0, 0.3, 0.47. Although not recorded here, we carried out studies for q = 0.2 and 
q = 0.4 and found the behavior quite similar to that reported. 

The outcome, as described in the introduction, is that the Humpherys-Zumbrun al- 
gorithm, implemented as /u(x), polar, or polar radial method, outperforms an optimized, 
adaptive-mesh, version of the Lee-Stewart algorithm by factor 1-10, and outperforms the 
actual fixed-mesh version proposed in |LS] by factor 100-1000, even in the high activation 
energy /square- wave limit. The hybrid method did not perform significantly better in the 
square- wave limit, and performed significantly worse in other regimes, hence we recommend 
that this option be discarded. Likewise, the Evans method, valid only when c/>(2) << max</>, 
does not seem to perform as well as the basic fJ>(x) method (see Table [6]), and so does not 
seem to be worth the trouble of implementation in this special regime. However, this option 
should perhaps not be discarded without more systematic investigation than was carried 
out here. 
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(0.0001,4,10,0.0086,1) 


(0.0001, 40, 10,0. 018, lc+01) 
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Table 5: Table of computational statistics, (M,i?,points, relative error,time) , where — M is the value of nu- 
merical negative infinity, R is the radius of the contour used to compute, points is how many points 
the adaptive solver used to compute the mesh, err is the maximum relative error between two points 
of the output, and time is how long it took to compute the contour. Here 4>{u) = Ce £ ^ T '"' where 



T{u) = 1- (1.5 



and C is chosen so that z passes through 0.5 at x = — 7. Here we are using 



the n(x) method. We take R to be the maximum of 2£ and the best root fit with e E o+Ei/\+E 2 /\ 
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method of Lee and Stewart (adaptive mesh in x) 


37.0 


homogenous method of Lee and Stewart (adaptive mesh in x) 


35.8 


adjoint n(x) method 


14.8 


hybrid method 


27.6 


polar adjoint method 


13.5 


inhomogeneous Erpenbeck method 


40.1 


polar adjoint radial method 


20.8 


adjoint (i method 


18.2 


polar /Drury method 


13.5 


fj,(x) method 


14.5 


polar radial method 


20.9 


(i method 


17.9 


Evans function method 


21.2 


homogenous Erpenbeck method 


32.4 


centered inhomogeneous Erpenbeck method 


51.8 



Table 6: Computational time comparison on a fixed A mesh of 55 points with radius 10 and inner radius 10 4 
for various methods with E = 10, q = 0.3, <j>(u) = e~ £/T(ll) , T(u) = 1 - (1.5 - uf , C = 10 21£/40 
and relative and absolute tolerance is set at 10 -12 in the ODE solver. 
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Table 7: Efficiency comparison for several methods: A- the p.(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and Stewart 
using a 4th order fixed mesh Rungc Kutta solver, F- the polar adjoint method, and G- the polar adjoint 
radial method. Here £ = 0, q = 0, C = 10 21£ / 40 , and 4>(u) = Ce- £ / T ( u ), T(u) = 1 - (u - 1.5) 2 . 
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Table 8: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and Stewart 
using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar adjoint 
radial method. Here £ = 0, q = 0.3, C = 10 21£ / 40 , and </>(u) = Ce- £ ' T( ~ U \ T(u) = !-(«- 1.5) 2 . 
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Table 9: Efficiency comparison for several methods: A- the fi(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and Stewart 
using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar adjoint 
radial method. Here £ = 0, q = 0.47, C = 10 21£ / 40 , and cf>(u) = Ce- £ / T< - u \ T{u) = 1 - (u - 1.5) 2 . 





A = 1 






A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


A 


























B 


1.0c-2 


0.0027 


2 


l.Oc-2 


0.0028 


2 


l.Oc-2 


0.0026 


2 


1.0e-2 


0.0028 


2 


C 


1.0c-2 


0.01 


10 


l.Oc-2 


0.01 


9 


l.Oc-2 


0.021 


23 


1.0e-2 


0.037 


39 


D 


1.0c-8 


0.39 


485 


l.Oc-8 


0.39 


438 


1.0e-9 


4.9 


6134 


1.0e-9 


6.1 


6816 


E 


2.7c-6 


187 


524288 


2.7e-6 


210 


524288 






>lc+6 






>le+6 


F 


l.Oc-02 


0.0026 


2 


l.Oc-02 


0.0026 


2 


l.Oc-02 


0.0025 


2 


l.Oc-02 


0.0026 


2 


G 


l.Oc-02 


0.0040 


3 


l.Oc-02 


0.0041 


3 


l.Oc-02 


0.0045 


3 


l.Oc-02 


0.0050 


3 



Table 10: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint fJ.(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 10, q = 0, C = 10 21£ / 40 , and 0(«) = CeT £ l T ^, T(u) = 1 - (u - 1.5) 2 . 
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Table 11: Efficiency comparison for several methods: A- the fi(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and Stewart 
using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar adjoint 
radial method. Here £ = 10, q = 0.3, C = 10 21£ / 40 , and 4>{u) = Ce- £ / T( - u \ T(u) = 1 - (it - 1.5) 2 . 
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Table 12: Efficiency comparison for several methods: A- the fi(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and Stewart 
using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar adjoint 
radial method. Here £ = 10, q = 0.47, C = 10 21£ / 40 , and 4>{u) = Ce- £ / T ("), T(u) = 1 - (u - 1.5) 2 . 
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Table 13: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint /J.(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 40, q = 0, C = 10 21£ / 40 , and <f>(u) = Cer £ l T< ^, T(u) = 1 - (tt - 1.5) 2 . 
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Table 14: Efficiency comparison for several methods: A- the fi(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and Stewart 
using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar adjoint 
radial method. Here £ = 40, q = 0.3, C = 10 21£ / 40 , and <fi(u) = C e - f / T ( u ), T{u) = 1 - (tt - 1.5) 2 . 
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Table 15: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint /J.(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and Stewart 
using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar adjoint 
radial method. Here £ = 40, q = 0.47, C = 10 21£ / 40 , and 4>{u) = Ce- £ I T( - U \ T(u) = 1 - (it - 1.5) 2 . 
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0.052 


128 


1.8e-07 


1.5 


4096 


1.8e-07 


1.6 


4096 


F 


1.0e-03 


0.0092 


6 


l.Oc-05 


0.011 


9 


1.0e-02 


0.02 


21 


1.0e-04 


0.039 


38 


G 


1.0e-04 


0.0183 


18 


1.0e-04 


0.0206 


18 


1.0e-04 


0.067 


70 





Table 16: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 1, q = 0.1, C = e £ / 2 , and <f>(u) = Ce~ £ l u . 





A = 1 






A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


A 


1.0e-02 


0.011 


9 


1.0e-03 


0.013 


10 


1.0e-02 


0.028 


28 


1.0e-03 


0.044 


44 


B 


1.0e-04 


0.0099 


10 


1.0c-04 


0.011 


10 


1.0e-02 


0.024 


26 


1.0e-05 


0.086 


91 


C 


1.0e-02 


0.012 


10 


1.0e-04 


0.016 


13 


1.0e-02 


0.027 


30 


1.0e-04 


0.074 


79 


D 


1.0e-06 


0.036 


44 


l.Oe-08 


0.059 


66 


l.Oe-08 


0.43 


532 


l.Oe-08 


0.47 


532 


E 


6.4e-07 


0.092 


256 


7c-07 


0.1 


256 


7.6c-08 


2.9 


8192 


7.8c-08 


3.3 


8192 


F 


1.0e-04 


0.012 


11 


1.0e-04 


0.011 


10 


1.0e-04 


0.028 


32 


1.0e-05 


0.086 


91 


G 


1.0e-04 


0.0251 


25 


1.0e-04 


0.0280 


25 


l.Oc-05 


0.1018 


106 


l.Oe-05 


0.1326 


122 



Table 17: Efficiency comparison for several methods: A- the fi(x) method, B- the adjoint fj,(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 1, q = 0.3, C = e £ / 2 , and <f>(u) = Ce~ £ / U . 





A = 1 






A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


1 To1 


Time 


Mesh 


A 


1.0e-03 


0.014 


14 


1.0c-04 


0.021 


19 


1.0c-02 


0.035 


38 


1.0c-03 


0.075 


79 


B 


1.0e-05 


0.02 


23 


l.Oc-05 


0.024 


25 


1.0c-03 


0.034 


40 


l.Oc-05 


0.17 


185 


C 


1.0c-03 


0.016 


16 


l.Oc-05 


0.028 


28 


1.0c-02 


0.034 


40 


1.0c-04 


0.13 


137 


D 


1.0e-07 


0.063 


76 


l.Oe-08 


0.086 


96 


l.Oe-08 


0.67 


840 


l.Oe-08 


0.75 


854 


E 


1.5e-07 


0.37 


1024 


1.7c-07 


0.41 


1024 


2.5c-07 


5.9 


16384 


2.6c-07 


6.5 


16384 


F 


l.Oc-05 


0.02 


23 


1.0c-06 


0.035 


36 


1.0c-04 


0.044 


52 


l.Oc-05 


0.17 


185 


G 


1.0e-05 


0.0493 


51 


1.0c-06 


0.0783 


72 


1.0c-06 


0.2847 


297 


l.Oc-05 


0.3128 


290 



Table 18: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint fJ.(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 1, q = 0.47, C = e £ / 2 , and <j>(u) = Ce- £ / u . 





A = 1 






A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


A 


1.0c-02 


0.01 


8 


1.0e-04 


0.016 


13 


1.0e-02 


0.027 


27 


1.0e-03 


0.046 


45 


B 


1.0c-03 


0.0089 


8 


1.0e-04 


0.012 


10 


1.0e-03 


0.024 


26 


1.0c-04 


0.052 


55 


C 


1.0c-02 


0.0096 


8 


1.0c-04 


0.017 


14 


1.0c-02 


0.027 


29 


1.0e-04 


0.069 


74 


D 


1.0c-06 


0.038 


46 


l.Oc-07 


0.044 


49 


l.Oe-08 


0.4 


497 


l.Oe-08 


0.44 


500 


E 


3.3c-07 


0.092 


256 


3.3c-07 


0.1 


256 


5.6C-07 


1.5 


4096 


5.9c-07 


1.6 


4096 


F 


1.0c-03 


0.0092 


8 


1.0c-04 


0.012 


10 


1.0c-02 


0.024 


25 


1.0e-04 


0.053 


55 


G 


lOc-04 


0.0202 


20 


1.0e-04 


0.0240 


21 


1.0c-04 


0.0811 


84 


1.0c-04 


01229 


113 



Table 19: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint fJ.(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 10, q = 0.1, C = e £ l 2 , and <f>(u) = Ce~ £ l u . 
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A — 1 






A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


A 


l.Oe-03 


0.014 


13 


l.Oc-05 


0.035 


32 


l.Oc-02 


0.047 


50 


l.Oe-04 


0.14 


155 


B 


1.0e-02 


0.0082 


8 


l.Oc-05 


0.027 


28 


l.Oc-03 


0.043 


49 


l.Oc-05 


0.21 


228 


C 


1.0e-02 


0.013 


12 


1.0c-05 


0.033 


32 


1.0c-02 


0.045 


52 


l.Oc-05 


0.24 


260 


D 


1.0e-07 


0.076 


93 


l.Oc-08 


0.1 


116 


l.Oc-08 


0.81 


1027 


l.Oc-08 


0.9 


1037 


E 


2.6e-07 


0.37 


1024 


3c-07 


0.41 


1024 


4.5C-07 




16384 


4.8c-07 


6.6 


16384 


F 


l.Oe-04 


0.013 


14 


l.Oc-06 


0.041 


42 


l.Oc-05 


0.062 


71 


1.0e-05 


0.21 


228 


G 


l.Oe-03 


0.031 


32 


l.Oc-05 


0.0528 


48 


1.0e-06 


0.2576 


269 


l.Oc-05 


0.3669 


339 



Table 20: Efficiency comparison for several methods: A- the fi(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 10, q = 0.3, C = e £ / 2 , and 4>(u) = Ce~ E l u . 





A = 1 








A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 




Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


A 


l.Oe-03 


0.029 


31 




l.Oc-08 


0.5 


540 


1.0e-02 


0.14 


174 


l.Oc-09 


7.3 


8009 


B 


1.0c-02 


0.019 


22 




l.Oc-08 


0.51 


542 


l.Oe-04 


0.16 


186 


l.Oc-08 


4.7 


5160 


C 


l.Oe-03 


0.028 


31 




l.Oc-08 


0.5 


542 


1.0e-02 


0.14 


176 


l.Oc-08 


4.7 


5169 


D 


l.Oc-08 


0.38 


482 




l.Oc-08 


0.42 


481 


l.Oc-09 


5.9 


7352 


l.Oc-09 


6.6 


7464 


E 


1.6e-07 


5.9 


16384 




1.4c-07 


6.6 


16384 














F 


l.Oe-03 


0.026 


29 




1.0e-06 


0.19 


201 


l.Oe-04 


0.21 


245 


1.0c-06 


1.8 


1933 


G 


1.0e-07 


0.3950 


411 




l.Oe-06 


0.6063 


565 


1.0e-07 


3.1675 


3218 


1.0c-06 


6.0556 


5575 



Table 21: Efficiency comparison for several methods: A- the fi(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 10, q = 0.47, C = e £ / 2 , and <f>(u) = Ce~ e / U . 





A = 1 






A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


A 


1.0c-02 


0.011 


9 


l.Oc-05 


0.022 


19 


1.0e-02 


0.032 


33 


l.Oe-04 


0.08 


84 


B 


1.0e-02 


0.0081 


7 


l.Oc-05 


0.016 


16 


l.Oe-04 


0.031 


35 


l.Oc-05 


0.11 


113 


C 


1.0e-02 


0.011 


10 


l.Oc-05 


0.021 


20 


l.Oc-02 


0.033 


36 


l.Oc-05 


0.14 


148 


D 


1.0e-06 


0.041 


50 


1.0c-07 


0.048 


54 


l.Oc-08 


0.48 


595 


l.Oc-08 


0.53 


597 


E 


8.1e-08 


0.18 


512 


8c-08 


0.21 


512 














F 


l.Oe-03 


0.0083 


8 


l.Oc-03 


0.0092 


8 


l.Oc-02 


0.027 


30 


l.Oe-04 


0.07 


74 


G 


l.Oc-04 


0.0232 


23 


l.Oe-04 


0.0282 


25 


l.Oc-05 


0.1116 


115 


l.Oc-05 


0.1686 


155 



Table 22: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint fJ.(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 20, q = 0.1, C = e £ / 2 , and <f>{u) = Ce~ £ / U . 





A = 1 






A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


A 


l.Oe-03 


0.021 


22 


l.Oc-09 


0.4 


418 


l.Oc-02 


0.093 


106 


l.Oc-09 


3.5 


3852 


B 


l.Oe-04 


0.022 


26 


l.Oc-09 


0.41 


424 


l.Oe-03 


0.091 


107 


l.Oc-08 


2.3 


2529 


C 


1.0c-04 


0.026 


30 


l.Oc-09 


0.4 


423 


l.Oc-02 


0.09 


108 


l.Oc-08 


2.3 


2532 


D 


1.0e-08 


0.21 


265 


l.Oc-08 


0.22 


246 


l.Oc-09 


2.9 


3668 


l.Oc-09 


3.2 


3705 


E 


1.8e-07 


1.5 


4096 


1.6C-07 


1.6 


4096 














F 


l.Oe-03 


0.017 


19 


l.Oc-06 


0.094 


100 


l.Oe-04 


0.11 


128 


l.Oc-05 


0.54 


583 


G 


l.Oe-05 


0.0874 


90 


1.0c-06 


0.1626 


149 


l.Oc-06 


0.7249 


740 


l.Oc-06 


1.3564 


1255 



Table 23: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint fJ.(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here £ = 20, q = 0.3, C = e £ / 2 , and (j>(u) = Ce- £ / u . 
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A — 1 






A = i 






A = 10 






A = lOi 








Tol 


Time 


pnts 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


Tol 


Time 


Mesh 


A 


1.0e-03 


0.13 


156 


1.0e-09 


7.8 


8654 


1.0e-02 


1.2 


1412 


l.Oe-10 


1.3e+02 


136649 


B 


1.0e-05 


0.18 


221 


1.0c-09 


8 


8670 


1.0e-04 


1.2 


1430 


l.Oe-10 


1.3e+()2 


136946 


C 


1.0e-04 


0.15 


175 


1.0e-09 


7.9 


8657 


1.0e-03 


1.2 


1419 


l.Oe-10 


1.3e+02 


136948 


D 








l.Oe-10 


319.6 


258949 








1.0e-ll 


9270.4 


3.9e06 


E 


B.9e-07 


96 


262144 


5.2o-08 


l.le+02 


262144 














F 


1.0e-04 


0.18 


218 


1.0e-07 


2.7 


2882 


1.0e-03 


2.1 


2023 


1.0e-07 


27 


28684 


G 


1.0e-07 


5.8414 


5911 


1.0e-07 


18.2622 


16857 















Table 24: Efficiency comparison for several methods: A- the fJ.(x) method, B- the adjoint fi(x) method, C- the hybrid 
method, D- the method of Lee and Stewart (using an adaptive mesh in z), E- the method of Lee and 
Stewart using a 4th order fixed mesh Runge Kutta solver, F- the polar adjoint method, and G- the polar 
adjoint radial method. Here Z = 20, q = OAT, C = e £ / 2 , and 4>{u) = Ce~ £ / U . 
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